Structure of POPC Lipid Bilayers in OPLS3e Force Field

It is crucial for molecular dynamics simulations of biomembranes that the force field parameters give a realistic model of the membrane behavior. In this study, we examined the OPLS3e force field for the carbon–hydrogen order parameters SCH of POPC (1-palmitoyl-2-oleoylphosphatidylcholine) lipid bilayers at varying hydration conditions and ion concentrations. The results show that OPLS3e behaves similarly to the CHARMM36 force field and relatively accurately follows the experimentally measured SCH for the lipid headgroup, the glycerol backbone, and the acyl tails. Thus, OPLS3e is a good choice for POPC bilayer simulations under many biologically relevant conditions. The exception are systems with an abundancy of ions, as similarly to most other force fields OPLS3e strongly overestimates the membrane-binding of cations, especially Ca2+. This leads to undesirable positive charge of the membrane surface and drastically lowers the concentration of Ca2+ in the surrounding solvent, which might cause issues in systems sensitive to correct charge distribution profiles across the membrane.


Impact of TIP3P water model in OPLS3e
To clarify if water model has an impact on cation binding affinity of membranes in OPLS3e, we performed simulations using the TIP3P water model instead of SPC with 1000 mM concentration of NaCl or CaCl2, and without ions. Without ions (Fig. S1), and at 1000 mM concentrations (Fig. S2), TIP3P produced similar order parameters as SPC, suggesting that there is no major difference between SPC and TIP3P water models in OPLS3e with or without ions. Figure S1. Order parameters at full hydration for headgroup, backbone, and acyl chains in OPLS3e with SPC and TIP3P water models. Experimental values for the POPC ( 1 H-13 C NMR) at 300 K are from Ref. 1

Equilibration times of the ion binding and impact on order parameters
Our simulations show that quite long equilibration time is needed for ion binding to the membrane in OPLS3e (both Ca 2+ and Cl -) as the number of ions keeps rising in the vicinity of membrane (Fig. S4) Figure S6. Change of order parameters in the POPC headgroup α and β segments in response to additional NaCl in OPLS3e force field with scaled and unscaled ions. Scaling was done with the factor of 0.75 for the charge of both Na + and Clions. Experimental values for DPPC ( 2 H NMR) at 323 K are from Ref. 3. The average C-H bond order parameters of R and S hydrogens were used to calculate the difference between the baseline and different salt concentrations.

Scaling the ion charge in OPLS3e
Over-binding of ions to membranes is a problem in most of the current force fields and based on our studies this is the case also with OPLS3e force field. One reason behind the over binding might be too high charge of the ions. In previous studies, it has been proposed that scaling the ion charge might help to overcome mistakes resulting from the lacking electronic polarization (6,7). In electronic continuum correlation (ECC) theory, based on quantum mechanical calculations, scaling factor of ions is 0.75. This scaling of ion charges has been able to improve monovalent ion binding in some force fields, but has not been sufficient for the divalent CaCl2.
To further investigate behavior of ions in OPLS3e, we scaled the charge of monovalent NaCl and divalent CaCl2 by factor 0.75, as in ECC theory, and performed simulations using 1000 mM concentrations with scaled charges. Scaling enhanced order parameters for the alpha and beta carbons in the presence of additional 1000 mM NaCl (Fig. S6) but was not sufficient for CaCl2 (Fig. S7). However, with scaled charged of CaCl2, ion binding to membrane equilibrates much faster than with unscaled charges (Fig. S8).  include also partial charges of the polar region of POPC phospholipids to build the ECC-POPC model based on the Lipid14 force field (they also did this for CHARMM36 with same scaling) (7). The resulting ECC-POPC model produces experimental order parameter responses to NaCl and even for divalent CaCl2 unlike any other current force field, and can thus be considered as one of the most realistic models to describe POPC lipid membrane in the presence of ions so far.

Small bug in CHARMM36 parameters
We noticed a small bug in the CHARMM36 force field parameters obtained at the time of the research (summer 2020) from CHARMM-GUI. Parameters were missing following dihedral line: 316a317 > CTL2 CEL1 CEL1 HEL1 9 1.800000e+02 2.510400e+01 2 By the end of 2020, the parameters in CHARMM-GUI were fixed. We compared the order parameters from simulations performed with parameters without (bug) and with the dihedral (bugfix) and discovered that there was no difference in order parameters between those parameter sets ( Fig S9). However, we have used CHARMM36 bugfix parameters in all the simulations presented in this paper.

Effect of temperature change in CHARMM36
In addition to 300 K, we simulated some CHARMM36 systems also using 314 K temperature. We can see that order parameters at full hydration are in general slightly higher at 314 K than 300 K (Fig. S10). Lowering hydration induces similar responses in both temperatures, but 314 K does not induce forking in 5 w/l hydration in contrast to 300 K (Fig. S11). Slight changes in order parameters might be linked to temperature-induced changes in area per lipid. At 314 K the area per lipid is larger, which allows lipids to move more than at 300 K (Fig. S12).

Area per lipid as a function of time
Area per lipid was used in this study to determine when the simulation systems had reached equilibrium. Other structural properties in addition to area per lipid (e.g., order parameters and density distributions) do not usually change significantly after area per lipid reaches a stable value (9). For simulations with ions we also waited for the ion distribution between the solvent and the membrane to stabilize which generally took longer than the stabilization of the area per lipid. In the case of CaCl2, we could not reach the equilibrium in the ion distribution in the 1000 ns production run. In most simulations area per lipid reaches a stable value during first 100 ns of the simulation (Fig S13.) Exceptions are systems with low hydration (5 w/l), and OPLS3e with additional CaCl2 for which area per lipid equilibrates significantly slower.